#script to run model for photosynthesis data – similar to growth_rate script

library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(afex)
## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: emmeans() for ANOVA models now uses model = 'multivariate' as default.
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
## 
## Attaching package: 'afex'
## The following object is masked from 'package:lme4':
## 
##     lmer
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(car)
library(MuMIn)
library (dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(emmeans)
library(DHARMa)
## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(performance)
library(patchwork)
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2
## ──
## ✔ tibble  3.1.8     ✔ purrr   1.0.1
## ✔ tidyr   1.3.0     ✔ stringr 1.5.0
## ✔ readr   2.1.3     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand()   masks Matrix::expand()
## ✖ rstatix::filter() masks dplyr::filter(), stats::filter()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ tidyr::pack()     masks Matrix::pack()
## ✖ dplyr::recode()   masks car::recode()
## ✖ purrr::some()     masks car::some()
## ✖ tidyr::unpack()   masks Matrix::unpack()
library(sjPlot)
## Learn more about sjPlot with 'browseVignettes("sjPlot")'.
library(sjmisc)
## 
## Attaching package: 'sjmisc'
## 
## The following object is masked from 'package:purrr':
## 
##     is_empty
## 
## The following object is masked from 'package:tidyr':
## 
##     replace_na
## 
## The following object is masked from 'package:tibble':
## 
##     add_case
library(mmtable2)
## 
## Attaching package: 'mmtable2'
## 
## The following object is masked from 'package:tidyr':
## 
##     table1
library(gt)
library(purrr)
library(stringr)
library(tidyr)

#load this file for rETR-based use per Silsbe and Kromkamp 2012 #all_runs_photosyn_data <- read.csv(“data_input/hyp_ulva_all_runs_ek_alpha.csv”)

#load this file for normalized to quantum efficiency of photosynthesis per same as above

all_runs_photosyn_data <- read.csv("data_input/hyp_ulva_all_runs_ek_alpha_normalized.csv")

assign run as a factor

all_runs_photosyn_data$Run <- as.factor(all_runs_photosyn_data$Run)

#assign temperature as a factor

all_runs_photosyn_data$Temperature <- as.factor(all_runs_photosyn_data$Temp...C.)

#assigns treatment as characters from integers then to factors

all_runs_photosyn_data$Treatment <- as.factor(as.character(all_runs_photosyn_data$Treatment))

assign deltaNPQ as a factor

all_runs_photosyn_data$deltaNPQ <- as.factor(all_runs_photosyn_data$deltaNPQ)

#toggle between the species for output. Use Day 9 for final analysis

ulva <- subset(all_runs_photosyn_data, Species == "ul" & RLC.Day == 9 & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol" 
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol" 
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol" 
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"

#ULVA pmax________________________________________________________________

ulva %>% ggplot(aes(pmax)) +
  geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

#run model without interaction between the treatments and temperature

ulva_pmax_model <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)

#construct null model to perform likelihood ratio test REML must be FALSE

ulva_pmax_treatment_null <- lmer(formula = pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
## boundary (singular) fit: see help('isSingular')
ulva_pmax_model2 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_pmax_treatment_null, ulva_pmax_model2)
## Data: ulva
## Models:
## ulva_pmax_treatment_null: pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_pmax_model2: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                          npar    AIC    BIC   logLik deviance  Chisq Df
## ulva_pmax_treatment_null    7 2042.7 2067.1 -1014.37   2028.7          
## ulva_pmax_model2           11 1919.9 1958.2  -948.94   1897.9 130.86  4
##                          Pr(>Chisq)    
## ulva_pmax_treatment_null               
## ulva_pmax_model2          < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ulva_pmax_temperature_null <- lmer(formula = pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_pmax_model3 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_pmax_temperature_null, ulva_pmax_model3)
## Data: ulva
## Models:
## ulva_pmax_temperature_null: pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_pmax_model3: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                            npar    AIC    BIC  logLik deviance Chisq Df
## ulva_pmax_temperature_null    9 1919.7 1951.0 -950.83   1901.7         
## ulva_pmax_model3             11 1919.9 1958.2 -948.94   1897.9 3.778  2
##                            Pr(>Chisq)
## ulva_pmax_temperature_null           
## ulva_pmax_model3               0.1512

#make residual plots of the data for ulva

hist(resid(ulva_pmax_model))

plot(resid(ulva_pmax_model) ~ fitted(ulva_pmax_model))

qqnorm(resid(ulva_pmax_model))
qqline(resid(ulva_pmax_model))

#check the performance of the model

performance ::check_model(ulva_pmax_model)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(ulva_pmax_model)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.5918687 0.6885467
tab_model(ulva_pmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  pmax
Predictors Estimates std. Error CI Statistic p df
(Intercept) 34.05 3.31 27.52 – 40.57 10.28 <0.001 229.00
Treatment [1] 21.54 3.54 14.57 – 28.51 6.09 <0.001 229.00
Treatment [2] 22.52 3.54 15.55 – 29.49 6.36 <0.001 229.00
Treatment [3] 41.51 3.54 34.54 – 48.49 11.73 <0.001 229.00
Treatment [4] 43.48 3.54 36.50 – 50.45 12.28 <0.001 229.00
Temperature [27] -5.55 2.80 -11.06 – -0.04 -1.99 0.048 229.00
Temperature [30] -3.22 2.54 -8.22 – 1.78 -1.27 0.205 229.00
Random Effects
σ2 135.25
τ00 Plant.ID 26.61
τ00 Run 8.13
τ00 RLC.Order 7.25
ICC 0.24
N Run 7
N Plant.ID 95
N RLC.Order 6
Observations 240
Marginal R2 / Conditional R2 0.592 / 0.689
plot(allEffects(ulva_pmax_model))

ulva %>% ggplot(aes(treatment_graph, pmax)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
  labs(x="salinity/nitrate", y= "Day 9 Pmax (μmols electrons m-2 s-1)", title= "A", subtitle = "Ulva lactuca") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
  ylim(-1, 170) + stat_mean() + 
  scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
  geom_hline(yintercept=50, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

#summarize the means for pmax

ulva %>% group_by(Treatment) %>% summarise_at(vars(pmax), list(mean = mean))
## # A tibble: 5 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          31.1
## 2 1          52.3
## 3 2          53.3
## 4 3          72.3
## 5 4          74.3
ulva %>% group_by(Run) %>% summarise_at(vars(pmax), list(mean = mean))
## # A tibble: 7 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      66.3
## 2 2      58.9
## 3 3      64.7
## 4 3.5    67.5
## 5 4      60.9
## 6 6      34.2
## 7 6.5    28.1
ulva %>% group_by(Treatment, RLC.Day) %>% summarise_at(vars(pmax), list(mean = mean))
## # A tibble: 5 × 3
## # Groups:   Treatment [5]
##   Treatment RLC.Day  mean
##   <fct>       <int> <dbl>
## 1 0               9  31.1
## 2 1               9  52.3
## 3 2               9  53.3
## 4 3               9  72.3
## 5 4               9  74.3

#add growth rate from other dataset to this one and subset by species

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$lunar.phase <- as.factor(growth_rate$lunar.phase)

#make a new column for weight change (difference final from initial)

growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100
gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)
ulva$lunar.phase <- (gr_ulva$lunar.phase)

#plot a regression between the photosynthetic independent variables of interest and growth rate

ulva_growth_etr_graph <- ggplot(ulva, aes(x=pmax, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Ulva lactuca Pmax vs Growth Rate", x = "Pmax (μmols electrons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_etr_graph
## `geom_smooth()` using formula = 'y ~ x'

#_____________________________________________________________________________________________________________ #HYPNEA

#There was no D9 RLC for hm6-4 on 11/12/21 but had to remove hm6-4 from 10/9/21 below to match growth data

hypnea <- subset(all_runs_photosyn_data, Species == "hm" & RLC.Day == 9 & uid != "2021-10-09_hm6-4")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol" 
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol" 
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol" 
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"

#for Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) # and hm6-4 on 10/9/21 because it was white and also looked dead

gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)

#hm1-1 on 10/12/22 and hm1-2 on 4/29/22 causing issues of influential observations #hm1-2 for both pmax and Ek – leaving them in dataset because no good reason to believe not good data #make a histogram and residual plots of the data for hypnea

hist(hypnea$pmax, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)

hypnea %>% ggplot(aes(pmax)) +
  geom_histogram(binwidth=5, fill = "#7D0033", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

#run model for pmax

hyp_pmax_model <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea)
hist(resid(hyp_pmax_model))

plot(resid(hyp_pmax_model) ~ fitted(hyp_pmax_model))

qqnorm(resid(hyp_pmax_model))
qqline(resid(hyp_pmax_model))

#check the performance of the model

performance::check_model(hyp_pmax_model)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(hyp_pmax_model)
##            R2m       R2c
## [1,] 0.2732395 0.6374791
summary(hyp_pmax_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: hypnea
## 
## REML criterion at convergence: 2271.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4349 -0.4862 -0.0761  0.4488  3.7551 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept)  42.71    6.535  
##  Run       (Intercept)  83.40    9.132  
##  RLC.Order (Intercept)  14.42    3.797  
##  Residual              139.86   11.826  
## Number of obs: 286, groups:  Plant.ID, 143; Run, 8; RLC.Order, 6
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)     56.175      7.120   6.156   7.890 0.000193 ***
## Treatment1       2.421      8.135   5.388   0.298 0.777153    
## Treatment2      -1.295      8.135   5.388  -0.159 0.879349    
## Treatment2.5    14.733     11.520   4.722   1.279 0.260159    
## Treatment3       3.905      8.135   5.388   0.480 0.650025    
## Treatment4       7.689      8.143   5.411   0.944 0.385299    
## Temperature27  -17.253      3.116  32.273  -5.537 4.05e-06 ***
## Temperature30  -19.429      2.729  52.871  -7.119 2.94e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.786                                          
## Treatment2  -0.786  0.956                                   
## Treatmnt2.5 -0.555  0.486  0.486                            
## Treatment3  -0.786  0.956  0.956  0.486                     
## Treatment4  -0.785  0.955  0.955  0.485  0.955              
## Temperatr27 -0.206  0.001  0.001  0.000  0.001  0.000       
## Temperatr30 -0.196  0.000  0.000  0.000  0.000  0.003  0.467
tab_model(hyp_pmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  pmax
Predictors Estimates std. Error CI Statistic p df
(Intercept) 56.18 7.12 42.16 – 70.19 7.89 <0.001 274.00
Treatment [1] 2.42 8.13 -13.59 – 18.44 0.30 0.766 274.00
Treatment [2] -1.29 8.13 -17.31 – 14.72 -0.16 0.874 274.00
Treatment [2.5] 14.73 11.52 -7.95 – 37.41 1.28 0.202 274.00
Treatment [3] 3.91 8.13 -12.11 – 19.92 0.48 0.632 274.00
Treatment [4] 7.69 8.14 -8.34 – 23.72 0.94 0.346 274.00
Temperature [27] -17.25 3.12 -23.39 – -11.12 -5.54 <0.001 274.00
Temperature [30] -19.43 2.73 -24.80 – -14.06 -7.12 <0.001 274.00
Random Effects
σ2 139.86
τ00 Plant.ID 42.71
τ00 Run 83.40
τ00 RLC.Order 14.42
ICC 0.50
N Run 8
N Plant.ID 143
N RLC.Order 6
Observations 286
Marginal R2 / Conditional R2 0.273 / 0.637
plot(allEffects(hyp_pmax_model))

#construct null model to perform likelihood ratio test REML must be FALSE

hypnea_pmax_treatment_null <- lmer(formula = pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_pmax_model2 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_pmax_treatment_null, hypnea_pmax_model2)
## Data: hypnea
## Models:
## hypnea_pmax_treatment_null: pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_pmax_model2: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                            npar    AIC    BIC  logLik deviance  Chisq Df
## hypnea_pmax_treatment_null    7 2335.5 2361.1 -1160.8   2321.5          
## hypnea_pmax_model2           12 2329.8 2373.7 -1152.9   2305.8 15.755  5
##                            Pr(>Chisq)   
## hypnea_pmax_treatment_null              
## hypnea_pmax_model2           0.007581 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hypnea_pmax_temperature_null <- lmer(formula = pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_pmax_model3 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_pmax_temperature_null, hypnea_pmax_model3)
## Data: hypnea
## Models:
## hypnea_pmax_temperature_null: pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_pmax_model3: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                              npar    AIC    BIC  logLik deviance  Chisq Df
## hypnea_pmax_temperature_null   10 2365.9 2402.4 -1173.0   2345.9          
## hypnea_pmax_model3             12 2329.8 2373.7 -1152.9   2305.8 40.109  2
##                              Pr(>Chisq)    
## hypnea_pmax_temperature_null               
## hypnea_pmax_model3            1.952e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hypnea %>% ggplot(aes(treatment_graph, pmax)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
  labs(x="Treatment", y= "Day 9 Pmax (μmols electrons m-2 s-1)", title= "B", subtitle = "Hypnea musciformis") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
  ylim(-1, 170) + stat_mean() + 
  scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
  geom_hline(yintercept=50, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

#plot a regression between the photosynthetic independent variables of interest and growth rate

hypnea_growth_etr_graph <- ggplot(hypnea, aes(x=pmax, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Hypnea musciformis Pmax vs Growth Rate", x = "Pmax (μmols electrons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor(label.y = 150)
hypnea_growth_etr_graph
## `geom_smooth()` using formula = 'y ~ x'

#summarize the means

hypnea %>% group_by(Treatment) %>% summarise_at(vars(pmax), list(mean = mean))
## # A tibble: 6 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          43.9
## 2 1          47.5
## 3 2          43.8
## 4 2.5        58.7
## 5 3          49.0
## 6 4          53.1
hypnea %>% group_by(Treatment, RLC.Day) %>% summarise_at(vars(pmax), list(mean = mean))
## # A tibble: 6 × 3
## # Groups:   Treatment [6]
##   Treatment RLC.Day  mean
##   <fct>       <int> <dbl>
## 1 0               9  43.9
## 2 1               9  47.5
## 3 2               9  43.8
## 4 2.5             9  58.7
## 5 3               9  49.0
## 6 4               9  53.1

#plot a regression between pmax and ek

hypnea_etr_ek_graph <- ggplot(hypnea, aes(x=pmax, y=ek.1)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Hypnea musciformis pmax vs Ek", x = "Pmax (μmols electrons m-2 s-1)", 
       y = "Ek (μmols photons m-2 s-1)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
hypnea_etr_ek_graph
## `geom_smooth()` using formula = 'y ~ x'